Estimation of Reference Evapotranspiration in a Semi-Arid Region of Mexico

Reference evapotranspiration (ET0) is the first step in calculating crop irrigation demand, and numerous methods have been proposed to estimate this parameter. FAO-56 Penman–Monteith (PM) is the only standard method for defining and calculating ET0. However, it requires radiation, air temperature, atmospheric humidity, and wind speed data, limiting its application in regions where these data are unavailable; therefore, new alternatives are required. This study compared the accuracy of ET0 calculated with the Blaney–Criddle (BC) and Hargreaves–Samani (HS) methods versus PM using information from an automated weather station (AWS) and the NASA-POWER platform (NP) for different periods. The information collected corresponds to Module XII of the Lagunera Region Irrigation District 017, a semi-arid region in the North of Mexico. The HS method underestimated the reference evapotranspiration (ET0) by 5.5% compared to the PM method considering the total ET0 of the study period (26 February to 9 August 2021) and yielded the best fit in the different evaluation periods (daily, 5-day mean, and 5-day cumulative); the latter showed the best values of inferential parameters. The information about maximum and minimum temperatures from the NP platform was suitable for estimating ET0 using the HS equation. This data source is a suitable alternative, particularly in semi-arid regions with limited climatological data from weather stations.


Introduction
Evapotranspiration (ET) is the sum of transpiration through the plant canopy and evaporation from the soil, plant, and free surface water [1,2]. ET is the most significant component of the hydrological cycle [1,3], due to which its estimation is of common interest in climatological, hydrological, forestry, and agricultural studies [4]. This last area ET is a fundamental variable for calculating water requirements, making efficient use of water in crop production [5]. ET can be measured directly using weighing lysimeters or by measuring the net flux of water vapor between the surface and the surrounding atmosphere using micrometeorological methods [6], which depend on the energy balance of the canopy and include the energy balance of Bowen's relation, eddy covariance, and the use of scintillometers [7].
Crop evapotranspiration (ET C ) is a crucial aspect of the water balance in agricultural areas. To estimate it, the most accessible method is to estimate reference evapotranspiration (ET 0 ) and then pair it with crop and soil coefficients [8]. Reference evapotranspiration (ET 0 ) is the evapotranspiration rate of a hypothetical reference crop (grass or alfalfa) with a height of 0.12 m, a fixed surface resistance of 70 s m −1 , and an albedo of 0.23, homogeneous, well-watered, free from diseases and pests, growing vigorously, and providing complete Module XII covers an area of 14,276.7 hectares with an elevation range of 1102 to 1114 m [56]. The slope of the area is gentle, at around 0.06%. It is oriented from south to north, with the southern part being the highest, as shown in Figure 2.  Module XII covers an area of 14,276.7 hectares with an elevation range of 1102 to 1114 m [56]. The slope of the area is gentle, at around 0.06%. It is oriented from south to north, with the southern part being the highest, as shown in Figure 2. Based on data gathered from Series VII of INEGI (2018) [57], it is estimated that the study area is primarily used for irrigated agriculture, with 89.3% of its surface area dedicated to this use. Human settlements make up 8.7% of the area, while the remaining 2% is used for other purposes (Table 1). The meteorological information used was daily averages for the period between 26 February (Julian day 57) and 9 August (Julian day 221) 2021 (n = 165). In the Lagunera Based on data gathered from Series VII of INEGI (2018) [57], it is estimated that the study area is primarily used for irrigated agriculture, with 89.3% of its surface area dedicated to this use. Human settlements make up 8.7% of the area, while the remaining 2% is used for other purposes (Table 1). The meteorological information used was daily averages for the period between 26 February (Julian day 57) and 9 August (Julian day 221) 2021 (n = 165). In the Lagunera Region Irrigation District 017, the main crops of the spring-summer cycle are grown in this period, including forage corn.
In addition, the meteorological variables were downloaded from the NP climate website (National Aeronautics and Space Administration-Prediction of Worldwide Energy Resource; https://power.larc.nasa.gov, accessed on 5 October 2022). This website collects information from various sources: data recorded on-site, satellite data, wind probes, and assimilated data systems [27].
The NNP weather data are based on a single assimilation model named GMAO (Global Modeling and Assimilation Office), starting from the MERRA-2 (Modern Era Retrospective-Analysis for Research and Applications) reanalysis dataset and the GEOS (Goddard Earth Observation System) data processing system [58,59]. Solar radiation is derived from the GEWEX SRB (Global Energy and Water Exchanges Project Surface Radiation Budget) project [60,61].
The horizontal resolution of the NP meteorological data source corresponds to a 1 2 • × 5/8 • latitude/longitude grid, and the solar data sources come from a 1 • × 1 • latitude/longitude grid. The current version no longer reassigns data to a common grid; once the data are processed and filed, they are available through the NP service package. The meteorological data is derived from NASA's GMAO MERRA-2 and GEOS 5.12.4 FP-IT. The NP platform team processes GEOS data daily and combines them with MERRA-2 data, producing daily time series that yield low-latency products usually available in approximately two days (real-time). Energy flow data (solar irradiance, thermal IR, and cloud properties) derive from NASA's GEWEX SRB Release 4-Integrated Product (R4-IP) file and CERES SYNIdeg and FLASHFlux projects. These data are processed daily and added to the daily time series, issuing products after approximately 4 days, almost in real-time [62].
The main features of the NP system database are shown in Table 2. The AWS is situated near the center of Module XII ( Figure 2) and aligns with the center of the NP platform cell. This suggests that one cell encompasses the entire study area's surface. Delayed data availability Approximately two days for temperature, relative humidity, and wind speed, and five days for solar radiation data.

ET 0 Estimation with Empirical Equations
ET 0 was estimated through three empirical equations with different information requirements: an equation based on a complete physical model (PM); another on temperature and solar radiation (HS); and the last one on temperature, relative humidity, and wind speed (BC). ET 0 was estimated daily with the FAO-56 Penman-Monteith method using Equation (1) [22,63]; this method is useful for arid, temperate, and tropical zones [19,22].
where R n is the net radiation at the reference crop surface (MJ m −2 d −1 ); G is the soil heat flux density (MJ m −2 d −1 ); u 2 is the wind speed at 2 m height (m s −1 ); e s is the saturation vapor pressure (kPa); e a is the actual vapor pressure (kPa); e s − e a is the vapor pressure deficit (kPa); ∆ is the slope of the vapor saturation pressure curve (kPa • C −1 ); T is the mean daily air temperature at 2 m height ( • C); and γ is the psychrometric constant (kPa • C −1 ). For daily time intervals, G values are relatively small, and therefore, this term was not included [22].

Hargreaves-Samani Method (ET 0 -HS)
The Hargreaves-Samani method estimates ET 0 based on temperature data only (Equation (2)) [64]. This equation was developed for semi-arid zones and is useful when solar radiation data are not available; however, as it is based on a few variables, its accuracy should be evaluated at the regional and local levels [65].
where T is the mean daily air temperature ( • C); R 0 is the extraterrestrial solar radiation (from tables, mm d −1 ); T max is the maximum daily air temperature ( • C); T min is the minimum daily air temperature ( • C); K H and K T are the empirical calibration parameters; and A H is a Hargreaves' empirical exponent. This study used the original values proposed by Hargreaves and Samani [64]: K H = 0.0023, K T = 17.78, and A H = 0.5.
where a and b are climatic calibration coefficients calculated with Equations (4) and (5), respectively; p is the mean annual percentage of daytime hours (value from tables, decimal); and T is the mean air temperature at 2 m height ( • C).
where RH MI N is the minimum relative humidity (%); n N is the ratio between theoretical and actual sunlit hours (value from tables, decimal).
where u 2 is the mean daily wind speed at 2 m height (m s −1 ). Table 3 shows the inferential parameters used for evaluating the empirical equations that estimate ET 0 (HS and BC), considering the PM method as a reference. Likewise, the climatic information of the NP platform was evaluated when calculating ET 0 through the reference method (PM NP ).

Inferential Evaluation Parameters
In the above equations, E i is the estimated value using the empirical equation; O i is the value obtained with the reference method; E is the average of estimated values obtained with the empirical equation; O is the average of the values obtained with the reference method; and n is the number of observations. The criteria for interpreting the reliability coefficient are cited in [67].
ET 0 was estimated in three different ways with empirical equations (daily, mean, and cumulative) using a total of 165 observations. The mean ET 0 was determined over the 5-day period, as was the cumulative. Table 3. Equations and optimal values of inferential parameters.

Parameter Equation Optimal
Value

Results and Discussion
The daily ET 0 calculated by the PM method and with AWS meteorological data ( Figure 3) had the peak value (8.8 mm d −1 ) on Julian day 126 (6 May 2021)-on the same day, a wind speed of 5.0 m s −1 was recorded, which was higher than the average recorded over the study period (2.2 m s −1 ). On the other hand, the minimum ET 0 (2.2 mm day −1 ) was recorded on Julian day 192 (11 July 2021)-the day that recorded a solar radiation value of 107.0 W m −2 , lower than the average for the study period (282.8 W m −2 ). This low radiation was due to atypical conditions: high cloudiness (rainfall of 7.6 mm recorded) and high relative humidity (84.5%). Some authors mention that wind speed and solar radiation are the climatic variables with the most significant influence on ET 0 estimates in the study area [27]. Other authors reach the same conclusion when performing a sensitivity analysis in other regions [68][69][70].
In the above equations, is the estimated value using the empirical equation; is the value obtained with the reference method; is the average of estimated values obtained with the empirical equation; is the average of the values obtained with the reference method; and is the number of observations. The criteria for interpreting the reliability coefficient are cited in [67].
ET0 was estimated in three different ways with empirical equations (daily, mean, and cumulative) using a total of 165 observations. The mean ET0 was determined over the 5day period, as was the cumulative.

Results and Discussion
The daily ET0 calculated by the PM method and with AWS meteorological data (Figure 3) had the peak value (8.8 mm d −1 ) on Julian day 126 (6 May 2021)-on the same day, a wind speed of 5.0 m s −1 was recorded, which was higher than the average recorded over the study period (2.2 m s −1 ). On the other hand, the minimum ET0 (2.2 mm day −1 ) was recorded on Julian day 192 (11 July 2021)-the day that recorded a solar radiation value of 107.0 W m −2 , lower than the average for the study period (282.8 W m −2 ). This low radiation was due to atypical conditions: high cloudiness (rainfall of 7.6 mm recorded) and high relative humidity (84.5%). Some authors mention that wind speed and solar radiation are the climatic variables with the most significant influence on ET0 estimates in the study area [27]. Other authors reach the same conclusion when performing a sensitivity analysis in other regions [68][69][70].   Table 4 shows the monthly and total ET 0 estimated using the empirical equations and the reference method (PM). Considering the month with the maximum ET 0 (May, with the HS and PM_NP equations and June with the BC method) and the reference method (PM), HS yielded an ET 0 value that was 6.6% lower vs. PM; BC, 12.5% lower; and PM_NP, 13.2% higher. However, considering the month with the minimum ET 0 (February) and the PM method, HS yielded an ET 0 12.8% higher vs. PM; BC, 6.8% higher; and PM_NP, 14.2% higher. It is observed that HS and BC underestimate ET 0 over most of the study period, consistent with the findings reported by some authors for an agroclimatic region similar to the study area [71]. However, when considering total ET 0 (whole study period) and the PM method, HS recorded an ET 0 value 5.5% lower vs. PM; BC, a value 15.6% lower; and PM_NP, 10.6% higher; therefore, HS was the equation that yielded values closest to the PM method. This is because HS considers temperature and radiation as the main energy sources that promote evapotranspiration [9,27].

Comparison of ET 0 Estimated by Empirical Equations versus the Reference Method
The results in Table 4 indicate an overestimation of ET 0 relative to the value obtained with the PM_NP method during the study period. The magnitude of this overestimation is related to the accuracy of each variable and has been reported only when using NP (NASA-POWER) data and the PM method [52,53,72]. Table 5 summarizes the relationship between the climatic variables recorded by the AWS and those obtained from the NP platform during the study period, where wind speed (WS) and solar radiation (SR) showed a low and moderate relationship, respectively. This same behavior has been reported by some authors for WS [27,45,73] and SR [74]. By contrast, Tmax and RH recorded a high ratio, and Tmin recorded a very high ratio. Some authors reported similar R 2 values for Tmin, Tmax [58], and RH [27] to those obtained in the present study. WS was the variable that yielded the lowest R 2 . This highlights the multiple challenges in determining this variable; these include quality control of the measured data since improving this aspect may return more accurate estimates [75]. Table 5. Relationship between the meteorological variables recorded by the automated weather station (AWS) and obtained from the NP platform during the study period.  Figure 4 depicts the bias in the data recorded by the automated weather station (AWS) relative to NP platform data for the following meteorological variables: temperature (maximum and minimum), relative humidity, solar radiation, and wind speed. It is observed that 44% of the maximum temperature data evaluated (n = 165) were virtually unbiased, while 39% of NP data overestimated Tmax by 2.1 • C to 7.5 • C, and the rest of the data (17%) underestimated Tmax by 1.2 • C to 5.5 • C (Figure 4a). Regarding the minimum temperature, 39% of the data evaluated showed bias values close to zero, while 46% of the NP data overestimated Tmin by 1.6 • C to 5.1 • C and the rest (15%) underestimated Tmin by 1.8 • C to 5.3 • C (Figure 4b).
Finally, 50% of the WS data showed bias values close to zero. It is also observe most data (41%) overestimated WS by 1.4 m s −1 to 2.9 m s −1 , and the rest of the data underestimated WS by 1.2 m s −1 to 3.3 m s −1 (Figure 4e).
Based on the above, the NP platform tends to overestimate Tmax, Tmin, SR, an while it underestimates RH. This same behavior was reported by Jiménez et al., study area for Tmin, WS, and RH [27].  The less biased RH values (values close to zero) were observed in 16% of the evaluated data; the NP platform underestimated RH by 3.0% to 38.8% in 76% of the data, and the rest of the data (7%) overestimated RH by 5.9% to 14.9% (Figure 4c). On the other hand, 45% of the evaluated data showed the minimum differences in solar radiation bias (values between 0 MJ m −2 d −1 and 1.5 MJ m −2 d −1 ), while 36% of the data overestimated radiation by 3.7 MJ m −2 d −1 to 17.2 MJ m −2 d −1 , and the rest of the data (19%) underestimated radiation by 2.9 MJ m −2 d −1 to 9.7 MJ m −2 d −1 (Figure 4d).
Finally, 50% of the WS data showed bias values close to zero. It is also observed that most data (41%) overestimated WS by 1.4 m s −1 to 2.9 m s −1 , and the rest of the data (9%) underestimated WS by 1.2 m s −1 to 3.3 m s −1 (Figure 4e).
Based on the above, the NP platform tends to overestimate Tmax, Tmin, SR, and WS while it underestimates RH. This same behavior was reported by Jiménez et al., in the study area for Tmin, WS, and RH [27].
Estimating the 5-day cumulative ET 0 improved the values of R 2 , r, and c relative to daily ET 0 and 5-day mean ET 0 . This behavior is consistent with the one reported by Jiménez et al. [27], who obtained better R 2 and RMSE values when estimating 10-day mean ET 0 versus daily data. Also, this way of estimating ET 0 yielded reliability coefficients (c) rated as "very good" for BC and PM_NP and "good" for HS. However, PM_NP showed the best R 2 , r, and c values, while HS yielded the best RMSE, PE, MBE, and b ( Table 6). The latter parameter returned values close to 1, indicating that the estimated values are statistically similar to observed or reference values [16]. Some authors reported similar RMSE values (1.1 mm d −1 ) when comparing ET 0 estimated by the HS equation and the PM method on a daily basis [71,76]. However, some authors recorded an RMSE (0.7 mm d −1 ) for 10-day mean data, which is similar to the RMSE value obtained in the present study for 5-day cumulative ET 0 [27]. Table 6. Comparison and inferential parameters to determine the ET 0 equation that best fits the study area. When graphically comparing the empirical equations versus the reference method (PM), daily ET 0 and 5-day mean ET 0 show a greater variability (Figure 5a,b); the 5-day cumulative ET 0 returned the best fit, with a lower variability of ET 0 values between the empirical equations and the PM method (Figure 5c).  Table 7 shows the results of the goodness-of-fit tests between ET0 calculated by the PM, HS, and BC methods, using maximum and minimum temperature data from the NP platform for the different calculation periods (daily, 5-day mean, and 5-day cumulative). The analyses of variance, with a 95% confidence interval (p-value < 0.0001), indicate a significant linear relationship between the PM method and the HS_NP and BC_NP equations for the three calculation periods. HS_NP yielded the best values of inferential parameters versus BC_NP, except for R 2 and r, in the daily ET0 estimate. This indicates that HS with NP In addition, HS yielded a better fit than the reference method (PM) in the three ways of estimating ET 0 . Other authors have also reported a better fit with the HS equation relative to other methods and have taken PM as a reference for arid and semi-arid regions [16,35,76]. This equation underestimates ET 0 over most of the study period because the methods based on solar temperature and radiation do not include wind speed [19]. Table 7 shows the results of the goodness-of-fit tests between ET 0 calculated by the PM, HS, and BC methods, using maximum and minimum temperature data from the NP platform for the different calculation periods (daily, 5-day mean, and 5-day cumulative). The analyses of variance, with a 95% confidence interval (p-value < 0.0001), indicate a significant linear relationship between the PM method and the HS _NP and BC _NP equations for the three calculation periods. HS _NP yielded the best values of inferential parameters versus BC _NP , except for R 2 and r, in the daily ET 0 estimate. This indicates that HS with NP temperature data is a suitable option for estimating ET 0 for different periods. In addition, we found that the estimation percent error (PE) is lower than 5% with HS _NP for the three ET 0 calculation periods. In addition, MBE is negative in the three periods, pointing to an underestimation with the HS _NP method. Some authors report this same ET 0 underestimation effect in semi-arid regions during the winter-summer period [27,71]. The 5-day mean ET 0 estimate recorded the best values in most statistical parameters relative to the mean daily ET 0 . However, the 5-day cumulative ET 0 estimate recorded the best r and R 2 values compared with the other two estimates (daily ET 0 and 5-day mean ET 0 ). These good results are obtained because grouping ET 0 over five days mitigates the variation in daily temperature associated with precipitation, wind speed, and cloudiness [77]. Figure 6 shows the dispersion of the calibrated HS method (HS _cal ) relative to PM for the different ET 0 calculation periods, depicting the best data fit obtained using data accumulated over five days.

Comparison of Estimated ET 0 with Observed (AWS) versus Estimated (NP) Data
The comparison of ET 0 estimates with the HS equation using temperature data recorded by the AWS (HS AWS ) and from the NP platform (HS NP ) yielded a high correlation for the three estimates ( Figure 7); the 5-day cumulative ET 0 recorded the highest R 2 . These R 2 values indicate the feasibility of estimating ET 0 using the NP platform's temperature data and the HS formula [27,71,76].
The accuracy of ET 0 calculated with methods BC and HS was compared versus the FAO-56 Penman-Monteith (PM) reference method using data from an automated weather station (AWS) and the NASA-POWER platform (NP). In this comparison, the HS equation returned the best fit in the different ways of estimating ET 0 : daily, 5-day mean, and 5-day cumulative, with the latter yielding the best fit.  Figure 6 shows the dispersion of the calibrated HS method (HS_cal) relative to PM for the different ET0 calculation periods, depicting the best data fit obtained using data accumulated over five days. The comparison of ET0 estimates with the HS equation using temperature data recorded by the AWS (HSAWS) and from the NP platform (HSNP) yielded a high correlation for  The accuracy of ET0 calculated with methods BC and HS was compared versus the FAO-56 Penman-Monteith (PM) reference method using data from an automated weather station (AWS) and the NASA-POWER platform (NP). In this comparison, the HS equation returned the best fit in the different ways of estimating ET0: daily, 5-day mean, and 5-day cumulative, with the latter yielding the best fit.

Conclusions
The Hargreaves-Samani (HS) method underestimated by 5.5% the reference evapotranspiration (ET 0 ) compared to the FAO-56 Penman-Monteith (PM) method considering the total ET 0 of the study period (26 February to 9 August 2021). This was because the calculation of ET 0 with the HS equation does not consider wind speed, which influences the evapotranspiration rate sometimes during the year in the study area. Nonetheless, this method is an alternative for calculating ET 0 in semi-arid regions for which only temperature records are available.
The HS equation yielded the best estimate relative to the reference method (PM) in the different ways of estimating ET 0 during the spring-summer crop cycle; the 5-day cumulative ET 0 showed the best fit. Therefore, this method is suitable for use with remotesensing data to determine crop evapotranspiration (ET c ) with 5-day temporal resolution images. It is necessary to conduct testing of HS in various agroclimatic conditions and perform a regional spatial evaluation using data from additional automated weather stations, such as those located within the entire 017 irrigation district.
The maximum and minimum temperature data from the NASA-POWER (NP) platform was suitable for estimating ET 0 with the HS equation. This data source is a timely alternative, particularly in semi-arid regions without data from weather stations.
The results showed that NP is a reliable data source for programming medium-and low-frequency irrigation (sprinkler and surface irrigation), which are common in the study area. In addition, they provide spatially comprehensive data, unlike the point values recorded by weather stations, which could be an enormous advantage when studying large regions, such as irrigation districts.  Data Availability Statement: The data are not publicly available because they are currently used in an ongoing thesis.